# Loading in Packages
library(tidyverse)
library(ggtext)

here::i_am("notebook/EEG_ROPE_Sims.Rmd")
here::here()
## [1] "D:/source/repo/psy-project"
# loading in User-Defined Functions
source(here::here("my-scripts/R", "eeg_rope_pipeline.R"))
source(here::here("my-scripts/R", "eeg_sim_functions.R"))

Chapter 4 of Project - Simulated EEG and ROPE Application

Overview of Chapter 1 - Multivarient Approach

Part One - Generate Raw EEG data

In this step, we use 1/F noise generated from a functions derived from Hank Steven’s work (original function is available on GitHub). After define to templates, one containing an effect and one without, we combine these with noise in order to represent a single trial. We then regenerate the noise again for each participant.

Defining Parameters to Create the Noise and Templates

set.seed(314)

.true_onset <- 450  ## Effect at 450ms
.stim_on <- 300     ## Stimulus is shown at 300ms
.freq <- 4         ## Defines resolution of stimulation, frequency of electrode (This is around 250 reading a second)
.max_time <- 800    ## End time for Epoch

.Xf <- seq(0, .max_time, .freq) ## start from zero and end at max time, going up in 4ms (freq)
Nf <- length(.Xf)   ## Number of timepoints

# Template one only contains zeros
.temp1 <- vector(mode = "numeric", length = Nf)  

# define magnitude of large effect
.magnit <- 2

# Generate ERP Peak
.erp <- dnorm(seq(-1.5, 1.5, length.out= 200/.freq), 0, 1)
.erp <- .erp - min(.erp)
.erp_small <- .erp / max(.erp)
.erp_large <- .magnit*.erp_small

# Getting length and storing them for later use
.l_erp <- length(.erp)
.l_pre_stim <- ceiling(.true_onset / .freq)

# Template 2 contains the ERP peak
temp1 <- c(rep(0, .l_pre_stim), .erp_small, rep(0, (Nf - .l_erp - .l_pre_stim)))
temp2 <- c(rep(0, .l_pre_stim), .erp_large, rep(0, (Nf - .l_erp - .l_pre_stim)))

Visualising the ERP Peak

tibble(x = .Xf, y = temp2) |> 
  ggplot(aes(x, y)) +
  geom_line()+
  geom_line(aes(y = temp1))+
  annotate("rect", xmin = .true_onset, xmax = (.true_onset+200)-.freq, ymin = 0, ymax = .magnit,
           alpha = .5, fill = "forestgreen")+
  annotate("rect", xmax = .true_onset, xmin = 0, ymin = 0, ymax = .magnit, alpha = .3)+
  annotate("rect", xmax = .max_time, xmin = (.true_onset+200)-.freq, ymin = 0, ymax = .magnit, alpha = .3)+
  annotate("text", label = "Null Hypothesis", x = 10, y = .magnit * .97, hjust = 0)+
  annotate("text", label = "Alternative Hypothesis", x = .true_onset, y = .magnit * .97, hjust = 0, colour = "forestgreen")+
  theme_minimal()+
  scale_y_continuous("Time (ms)", expand = c(0,0))

Combining Templates with Noise - Simulating 12 “experiments”, starting with 2 trials and doubling each time (until 4096 trials)

sample_size <-c(10, 25, 50, 100, 150, 200) # number of trials
.gsp <- 1 # gamma spectral power 
.outvar <- 1 # noise variance

alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE

effect_time <- c(.true_onset - .stim_on, .true_onset + 200 - .freq - .stim_on)
time_cond <- expression((time < effect_time[1] | time > effect_time[2]))

## storing values in list, so make functions neater
eeg_pipeline_attr <- list(
  sample_size = sample_size, # Defines Sample Size/ Number of Trials 
  alpha = alpha, 
  num_time_points = Nf, 
  max_time = .max_time, 
  cond1_base = temp1, 
  cond2_base = temp2, 
  gamma_spec_power = .gsp, 
  noise_var = .outvar,   # Influence Noise output - Zero means no Noise
  stim_onset = .stim_on  # Time when Stimulus was "shown"
)

## No effect Present (Only using temp1)

eeg_pipeline_attr2 <- eeg_pipeline_attr
eeg_pipeline_attr2$cond2_base <- eeg_pipeline_attr$cond1_base 
## Effect Present in template 2 (temp2)
sim_df_1 <- sim_eeg_pipeline(eeg_pipeline_attr, seed = seed)

sim_df_2 <- sim_eeg_pipeline(eeg_pipeline_attr2, seed = seed)
## df contains p, t, ci and rope min and max values, will also be used to plot line of significance at bottom of graphs
rope_df_1 <- sim_df_1 |> 
  generate_rope_df(eeg_pipeline_attr, 
                   rope_func = get_eeg_nu_rope,
                   static_margins = static_margins, 
                   replication = F, 
                   sequential = F, 
                   null_rope = T)
## 
## A Null ROPE was used to Define Significance
rope_df_2 <- sim_df_2 |> 
  generate_rope_df(eeg_pipeline_attr2, 
                   rope_func = get_eeg_nu_rope,
                   static_margins = static_margins, 
                   replication = F, 
                   sequential = F, 
                   null_rope = T)
## 
## A Null ROPE was used to Define Significance
.p1 <- sim_df_1 |> 
  plot_sim_base(rope_df = rope_df_1) 

.p2 <- sim_df_2 |> 
  plot_sim_base(rope_df = rope_df_2)

.p1
## [[1]]

.p2
## [[1]]

# ggsave(plot = p1, filename = here::here("images", "ch4_rope_eff_plain.png"), width = 10, height = 8, dpi = 360)
# ggsave(plot = p1, filename = here::here("images", "ch4_rope_eff_plain.png"), width = 10, height = 8, dpi = 360)

Part Two - Define Nulls ROPEs

Experiments with an Effect

We are going to define two ROPE (Regions of Practical Equivalence) using the data to represent the null and alternative hypothesis. The Null ROPE will be defined by the 90% quantile range of the mean pre-stimulus electrical activity for both groups and the Alternative ROPE will be defined by the mean differences of only the significant p-values. It may be better specify the range of time points a priori, however, I will use two different ROPEs to represent the alternative hypthesis, one defined from raw p-values and the other by bonferroni corrected p-values. This is provide a more conservative definition of an effect, as some results pre-stimulus can be significant and those false positive should not influence the alternative ROPE.

Once the ROPE is created for an “experiment”, it will superimpose the data of the next experiment, to see how the new data fits to the previous study. If the Confidence intervals are perfectly contained in a ROPE, then we can consider that as support for rejecting and maintaining the null hypothesis. For example, if the Confidence interval is inside the Null ROPE, even if the difference is statistically significant, this should likely be interpreted as evidence for keeping the null. Likewise, if the alternative ROPE covers the majority of the ERP amplitude, then that would be evidence of a replication of the underlying effect and thus, evidence to reject the null.

.ch4_sim_1 <- plot_rope(.p1, rope_df_1, null_rope = T, change_subtitle = F)

.ch4_sim_1
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

As sample size increases the effect becomes more clear and obvious, and the background noise appears to cancel itself out. This in turn increases the size of the alternative ROPEs, both ROPE defined from raw p-values and bonferroni corrected are wided at the larger sample size. However, raw p-values cause the ROPE to cover the entire range of values, including values where an effect is unlikely. Moreover, it appears only the Bonferroni defined ROPE maintains a gap between the two ROPEs, the Null and Alternative, whch would allow easier interpretation. For, if the ROPEs overlap, and the Confidence interval stretches across both, then it is unclear which ROPE it supports, preventing us to make a clear preference of a hypothesis.

It is worth noting that the peak of each of the effects appears to go above the ROPE, as the quantile range which is used to define these ROPEs would need to contain values larger than what is observed. It may not be suitable to use a ROPE to investigate the maximum values, instead it should likely be used to prevent false positives and ensure the observed difference is outside an interval which represents background noise.

Experiments without an Effect

To see how a ROPE can prevent false positives, we ran the EEG simulation but with no underlying effect, therefore any statistically significant results are false positives.

.ch4_sim_2 <- plot_rope(.p2, rope_df_2, null_rope = T)

.ch4_sim_2
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

Replication of ERP peaks

rep_df_1 <- sim_repli_eeg_pipeline(pipeline_attr = eeg_pipeline_attr, seed = seed, total_experiments = 2)
rep_df_2 <- sim_repli_eeg_pipeline(pipeline_attr = eeg_pipeline_attr2, seed = seed, total_experiments = 2)
rep_rope_1 <- rep_df_1 |>
  generate_rope_df(pipeline_attr = eeg_pipeline_attr, 
                   rope_func = get_eeg_nu_rope,
                   static_margins = static_margins, 
                   replication = T, 
                   sequential = F, 
                   null_rope = T)
## 
## A Null ROPE was used to Define Significance
rep_rope_2 <- rep_df_2 |>
  generate_rope_df(pipeline_attr = eeg_pipeline_attr2, 
                   rope_func = get_eeg_nu_rope,
                   static_margins = static_margins, 
                   replication = T, 
                   sequential = F, 
                   null_rope = T)
## 
## A Null ROPE was used to Define Significance
# ropes <- unique(rep_rope_1$rope)
.rep_p1 <- rep_df_1 |>
  plot_rep_base(rep_rope_1)

.rep_p2 <- rep_df_2 |>
  plot_rep_base(rep_rope_2)
.rep_p1
## [[1]]

.rep_p2
## [[1]]

.ch4_rep_1 <- plot_rope(.rep_p1, rep_rope_1, null_rope = T, replication = T)
.ch4_rep_1[[1]]

.ch4_rep_2 <- plot_rope(.rep_p2, rep_rope_2, null_rope = T, replication = T)
.ch4_rep_2
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

### Alternative ROPE

rope_df_1a <- sim_df_1 |> 
  generate_rope_df(eeg_pipeline_attr, 
                   rope_func = get_eeg_alter_rope, 
                   replication = F, 
                   sequential = F, 
                   null_rope = F,
                   replicate.effect = T)
## 
## An Alternative ROPE was used to Define Significance
rope_df_2a <- sim_df_2 |> 
  generate_rope_df(eeg_pipeline_attr2, 
                   rope_func = get_eeg_alter_rope,
                   replication = F, 
                   sequential = F, 
                   null_rope = F, 
                   replicate.effect = T)
## 
## An Alternative ROPE was used to Define Significance
## adding alternative rope to previous simulations
.ch4_sim_1a <- plot_rope(.ch4_sim_1, df = rope_df_1a, null_rope = F, change_subtitle = T, replication = F)
.ch4_sim_1a
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## adding alternative rope to previous simulations
.ch4_sim_2a <- plot_rope(.ch4_sim_2, df = rope_df_2a, null_rope = F, change_subtitle = T)
.ch4_sim_2a
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

Applying Alternative ROPE to replications

rep_rope_1a <- rep_df_1 |>
  generate_rope_df(pipeline_attr = eeg_pipeline_attr, 
                   rope_func = get_eeg_alter_rope,
                   replication = T, 
                   sequential = F, 
                   null_rope = F,
                   replicate.effect = T)
## 
## An Alternative ROPE was used to Define Significance
rep_rope_2a <- rep_df_2 |>
  generate_rope_df(pipeline_attr = eeg_pipeline_attr2, 
                   rope_func = get_eeg_alter_rope,
                   replication = T, 
                   sequential = F, 
                   null_rope = F, 
                   replicate.effect = T)
## 
## An Alternative ROPE was used to Define Significance
.ch4_rep_1a <- plot_rope(.ch4_rep_1, df = rep_rope_1a, null_rope = F, change_subtitle = T, replication = T)
.ch4_rep_1a
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

.ch4_rep_2a <- plot_rope(.ch4_rep_2, df = rep_rope_2a, null_rope = F, change_subtitle = T, replication = T)
.ch4_rep_2a
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

Dynamic Alpha Rule - min(alpha, beta)

library(pwr)
df_dyn <- sim_df_1 |> 
  filter(!eval(time_cond)) |> 
  group_by(n_trial) |> 
  summarise(
    mean_peak = mean(mn_diff_13), 
    d = abs(mean_peak)/ sqrt(0.5 * sd (sd_cond1) + sd(sd_cond3)), 
    half_d = 0.5 * d,
    .groups = "drop"
    ) |> 
  rowwise() |> 
  mutate(
    beta_2 = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = half_d)$power), 
    new_alpha_2 = min(0.05, beta_2), 
    beta = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = d)$power), 
    new_alpha = min(0.05, beta)
    )

dyn_stats <- sim_df_1 |> 
  pivot_for_stats(replication = F) |> 
  inner_join(df_dyn, by = "n_trial") |> 
  mutate(
    sig_rule1 = if_else(p < new_alpha, 0, NULL),
    sig_rule2 = if_else(p < new_alpha_2, 0, NULL)
  )

dyn_stats
## # A tibble: 2,412 × 17
##     time n_trial compari…¹      p     t  diff    ci_l  ci_u mean_…²     d half_d
##    <dbl>   <dbl> <chr>      <dbl> <dbl> <dbl>   <dbl> <dbl>   <dbl> <dbl>  <dbl>
##  1  -300      10 12        0.287   1.10 0.524 -0.481   1.53   0.374 0.617  0.308
##  2  -300      10 13        0.0419  2.20 0.895  0.0369  1.75   0.374 0.617  0.308
##  3  -296      10 12        0.230   1.24 0.596 -0.415   1.61   0.374 0.617  0.308
##  4  -296      10 13        0.0221  2.50 1.00   0.162   1.85   0.374 0.617  0.308
##  5  -292      10 12        0.228   1.25 0.547 -0.377   1.47   0.374 0.617  0.308
##  6  -292      10 13        0.0345  2.29 0.843  0.0684  1.62   0.374 0.617  0.308
##  7  -288      10 12        0.130   1.59 0.672 -0.219   1.56   0.374 0.617  0.308
##  8  -288      10 13        0.0169  2.64 0.913  0.185   1.64   0.374 0.617  0.308
##  9  -284      10 12        0.184   1.39 0.609 -0.320   1.54   0.374 0.617  0.308
## 10  -284      10 13        0.0246  2.45 0.889  0.128   1.65   0.374 0.617  0.308
## # … with 2,402 more rows, 6 more variables: beta_2 <dbl>, new_alpha_2 <dbl>,
## #   beta <dbl>, new_alpha <dbl>, sig_rule1 <dbl>, sig_rule2 <dbl>, and
## #   abbreviated variable names ¹​comparison, ²​mean_peak
.p1[[1]] + 
  geom_point(data = dyn_stats, aes(x = time, y = sig_rule1 - .9), colour = "darkgreen", inherit.aes = F, size = 1, shape = 15)+
  geom_point(data = dyn_stats, aes(x = time, y = sig_rule2 - .8), colour = "darkred", inherit.aes = F, size = 1, shape = 15)
## Warning: Removed 12714 rows containing missing values (`geom_point()`).
## Warning: Removed 2319 rows containing missing values (`geom_point()`).
## Warning: Removed 2275 rows containing missing values (`geom_point()`).

dyn_rep_1 <- rep_df_1 |> 
  filter(!eval(time_cond), experiment == "A") |> 
  group_by(n_trial) |> 
  summarise(
    mean_peak = mean(mn_diff_13), 
    d = abs(mean_peak)/ sqrt(0.5 * sd (sd_cond1) + sd(sd_cond3)), 
    half_d = 0.5 * d,
    .groups = "drop"
    ) |> 
  rowwise() |> 
  mutate(
    beta = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = d)$power), 
    beta_2 = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = half_d)$power), 
    new_alpha = min(0.05, beta),
    new_alpha_2 = min(0.05, beta_2)
    )

dyn_rep_stats <- rep_df_1 |> 
  pivot_for_stats(replication = T) |> 
  inner_join(dyn_rep_1, by = "n_trial") |> 
  mutate(
    sig_rule1 = if_else(p < new_alpha, 0, NULL),
    sig_rule2 = if_else(p < new_alpha_2, 0, NULL)
  )

## min difference required to be significant with new alpha rule
min_diff <- convert_t_to_raw(sim_df_1, alpha = dyn_rep_1$new_alpha, 
                             t_values = qt(dyn_rep_1$new_alpha, 100, lower.tail = F), null_rope = F) |> 
  unlist()

min_diff2 <- convert_t_to_raw(sim_df_1, alpha = dyn_rep_1$new_alpha_2, 
                              t_values = qt(dyn_rep_1$new_alpha_2, 100, lower.tail = F), null_rope = F) |> 
  unlist()

min_diff_df <- tibble(
  n_trial = unique(sim_df_1$n_trial),
  min_diff_max = min_diff,
  min_diff2_max = min_diff2,
  min_diff_min = -min_diff,
  min_diff2_min = -min_diff2,
       ) |> 
  mutate(across(where(is.double), ~ if_else(!(abs(.) == Inf), ., NULL)))
 
dyn_rep_stats <- inner_join(dyn_rep_stats, min_diff_df, by = "n_trial")

.rep_p1[[1]] + 
  geom_point(data = dyn_rep_stats, 
             aes(x = time, y = sig_rule1 - 1.25, fill = comparison), 
             inherit.aes = F, size = 1, shape = 22, stroke = NA )+
  geom_point(data = dyn_rep_stats, aes(x = time, y = sig_rule2 - 1.5, fill = comparison), 
             inherit.aes = F, size = 1, shape = 22, stroke = NA )+
  geom_ribbon(data = dyn_rep_stats, aes(x = time, ymin = min_diff_min, ymax = min_diff_max), inherit.aes = F, alpha = fill_alpha)+
  geom_ribbon(data = dyn_rep_stats, aes(x = time, ymin = min_diff2_min, ymax = min_diff2_max), inherit.aes = F, alpha = fill_alpha)

### Bayes Factors (Testing/ Not Final)

sim_bf_1 <- sim_eeg_pipeline(eeg_pipeline_attr, seed = seed, get_bf = T) |> 
  mutate(
    bf_sig = if_else(bf_13 > 3, 0, NULL),
    bf_null_sig = if_else((1/bf_13) > 3, 0, NULL),
    )
.p_bf1 <- sim_bf_1 |> 
  plot_sim_base()
## Argument 'rope_df' was empty - No lines of Significance were plotted
.p_bf1[[1]] + 
  geom_point(data = sim_bf_1, aes(x = time, y = bf_sig -1 ),shape = 22, fill = fill_colours[1], stroke = NA, inherit.aes = F)+
  geom_point(data = sim_bf_1, aes(x = time, y = bf_null_sig - 1),shape = 22, fill = fill_colours[2], stroke = NA, inherit.aes = F)
## Warning: Removed 1006 rows containing missing values (`geom_point()`).
## Warning: Removed 617 rows containing missing values (`geom_point()`).

### Bayesian Parametre Estimation and HDI (Testing/Not Final)

sim_df_1
## # A tibble: 1,206 × 20
##     time n_trial mn_co…¹ mn_co…² mn_co…³ sd_co…⁴ sd_co…⁵ sd_co…⁶ mn_di…⁷ mn_di…⁸
##    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1  -300      10  -0.184   0.340   0.710   1.00    1.13    0.807   0.524   0.895
##  2  -296      10  -0.326   0.270   0.679   0.920   1.20    0.874   0.596   1.00 
##  3  -292      10  -0.278   0.269   0.564   0.815   1.12    0.833   0.547   0.843
##  4  -288      10  -0.304   0.368   0.609   0.848   1.03    0.687   0.672   0.913
##  5  -284      10  -0.290   0.319   0.599   0.835   1.11    0.785   0.609   0.889
##  6  -280      10  -0.233   0.265   0.614   0.755   1.09    0.833   0.498   0.847
##  7  -276      10  -0.328   0.219   0.639   0.815   1.20    0.892   0.546   0.967
##  8  -272      10  -0.441   0.190   0.620   0.850   1.20    0.864   0.631   1.06 
##  9  -268      10  -0.538   0.201   0.574   0.681   1.12    0.851   0.739   1.11 
## 10  -264      10  -0.527   0.204   0.649   0.640   0.962   0.765   0.731   1.18 
## # … with 1,196 more rows, 10 more variables: p_12 <dbl>, p_13 <dbl>,
## #   t_12 <dbl>, t_13 <dbl>, ci_l_12 <dbl>, ci_u_12 <dbl>, ci_l_13 <dbl>,
## #   ci_u_13 <dbl>, bon_13_sig <dbl>, raw_13_sig <dbl>, and abbreviated variable
## #   names ¹​mn_cond1, ²​mn_cond2, ³​mn_cond3, ⁴​sd_cond1, ⁵​sd_cond2, ⁶​sd_cond3,
## #   ⁷​mn_diff_12, ⁸​mn_diff_13

Just like the ROPEs for when there is an underlying effect, the ROPE constructed from raw p-values appears to be less precise and helpful. Here we can see how False positves influence the ROPE whilst no results were significant after bonferroni correction.

Evaluating the False Positives with respect to no ROPE, and the different types of ROPEs

The use of a ROPE should reduce the number of false positives, as if the Confidence interval of a statistically significant point is not completely outside the Null ROPE then no conclusions can be made.

It is probably worth replicating this pipeline (the above process) several times to assess behaviour in the long run - as consistency is difficult to assess from a few trials. These following graphs are influenced by randomness, as the noise is randomly generated, so repeating this pipeline should reduce this unwanted influence.

seeds = 1:1000 # the seeds that were used in the 1000 repetitions of the Pipeline

## Reading the Pipeline Simulations - 1000 repetitions
sim_rope_1 <- read_csv(here::here("sim_data", "sim_egg_pipeline_eff.csv"))
sim_rope_2 <- read_csv(here::here("sim_data", "sim_egg_pipeline_no_eff.csv"))

## Evaluating the ROPEs of Number of False Positives
eval_rope_1 <- sim_rope_1 |> 
  eval_rope(alpha = 0.05, effect_time_vector = effect_time)

eval_rope_2 <- sim_rope_2 |> 
  eval_rope(alpha = 0.05, effect_time_vector = effect_time)
## Functions to create graphs
source(here::here("my-scripts", "R", "eeg_rope_graphs.R"))
.eval_rope_g1 <- eval_rope_1 |> 
  fp_graphs(effect = T)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

.eval_rope_g2 <- eval_rope_2 |> 
  fp_graphs(effect = F)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Raw p-values are going to generate false positives, which is expected after doing around 400 t-tests. The ROPE does seem to catch some of these, as if the Confidence interval overlaps with the ROPE, then we don’t conclude it is a significant difference and without a decision. However, bonferroni corrected p-values seems to be prefect as controlling for false positives, which is probably because of the multivariate approach used in the analysis. IF the EEG was analysed with the mean amplitude then bonferroni correction will start to resemble the raw p values, and hopefully the ROPE will still be effect.

Evaluating the True Positives with respect to no ROPE, and the different types of ROPEs

Whilst minimising False positives is desirable, we can achieve a false positive rate of zero by sacrificing our ability to detect meaningful statistical differences. This is the classic issue of Sensitivity vs Specificity. For example, if we set the margins of a ROPE to -10 and +10, we would have no False Positive because every value is inside that ROPE, so we would have no True positives.

.eval_rope_g3 <- eval_rope_1 |> 
  tp_graphs()
## [[1]]
## Warning: Removed 3 rows containing missing values (`geom_point()`).

## 
## [[2]]
## Warning: Removed 3 rows containing missing values (`geom_point()`).

## 
## [[3]]

## 
## [[4]]

# auc_df <- eval_rope_df |>
#   inner_join(tp_df_1, by = c("n_trial", "rope"))

## Not as neat as I expected - maybe we need more values or the data is not properly formatted?
# eval_rope_1|>
#   filter(n_trial != 10) |>
#   ggplot(aes(x = fpr, y = tpr, colour = rope))+
#   geom_line()+
#   theme_minimal()

Next steps: Repeat pipeline to get smooth curves and long-term behaviour of the various methods Repeat pipeline but analyse data using mean amplitude instead of t-tests at every time point Add references to Laken’s and Kruschke’s work, tidy up unprofessional language.

Overview of Chapter 2 - Mean Amplitude Approach

alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE

## defining region where peak occurs - used to find mean amplitude
time_cond <- expression((time < effect_time[1] | time > effect_time[2]))
## Process to proform anovas on each set of n_trial for both conditions
peak_anova_df_1 <- sim_df_1 |> 
  select(1:4) |>  # removing columns that are used for Multivarient Approach
  mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |> 
  pivot_longer(cols = starts_with("mn_cond"), names_to = "cond", values_to = "value") |> 
  mutate(cond = as.factor(cond)) |> 
  filter(!eval(time_cond)) |> 
  nest(data = -n_trial) |> 
  mutate(model = map(.x = data, ~ anova(lm(value ~ cond, .))), 
         tidy = map(model, broom::tidy)) |> 
  select(n_trial, tidy) |> 
  unnest(cols = c(tidy))
# 
# peak_df_1 <- sim_df_1 |> 
#   select(1:5) |>  # removing columns that are used for Multivarient Approach
#   mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |> 
#   group_by(n_trial, in_effect) |> 
#   summarise(
#     mn_peak = mean(mn_diff, na.rm = T),
#     mn_peak_t1 = mean(mn_diff, na.rm = T, trim = .2),
#     mn_peak_t1 = mean(mn_diff, na.rm = T, trim = .1),
#     mn_peak_t5 = mean(mn_diff, na.rm = T, trim = .05),
#     sd_peak = sd(mn_diff, na.rm = T), 
#     .groups = "drop"
#     ) |> 
#   pivot_longer(cols = contains("mn_"), 
#                names_to = "mean_type", 
#                values_to = "value")
  
peak_df_1 <- sim_df_1 |> 
  select(1:4) |>  
  mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |>
    dplyr::group_by(n_trial, in_effect) |> 
    dplyr::summarise(
      mn_peak_cond1 = mean(mn_cond1),
      mn_peak_cond2 = mean(mn_cond2),
      mn_diff = mn_peak_cond2 - mn_peak_cond1,
      p = t.test(mn_cond1, mn_cond2)[["p.value"]],
      ci_l = t.test(mn_cond2, mn_cond1, conf.level = 1 - alpha, alternative = "two.sided")[["conf.int"]][1],
      ci_u = t.test(mn_cond2, mn_cond1, conf.level = 1 - alpha, alternative = "two.sided")[["conf.int"]][2],
      .groups = "drop"
    ) |> 
    dplyr::group_by(n_trial) |>
    dplyr::mutate(raw_sig = dplyr::if_else(p < alpha, min(ci_l), NULL)) |>  # just so the graph has a line at zero to show significance
    dplyr::ungroup()



#   
# sim_df_2 |> 
#   select(1:5) # removing columns that are used for Multivarient Approach
#   
peak_df_1
## # A tibble: 12 × 9
##    n_trial in_eff…¹ mn_pea…² mn_pea…³  mn_diff       p     ci_l     ci_u raw_sig
##      <dbl>    <dbl>    <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>   <dbl>
##  1      10        0 -0.0991   0.0409   0.140   3.82e-4  0.0636   0.216    -0.594
##  2      10        1  0.899    0.465   -0.434   5.53e-7 -0.594   -0.274    -0.594
##  3      25        0 -0.0423   0.0156   0.0579  2.91e-3  0.0200   0.0958   -0.317
##  4      25        1  0.723    0.543   -0.180   1.09e-2 -0.317   -0.0423   -0.317
##  5      50        0 -0.00940 -0.0284  -0.0190  8.42e-2 -0.0406   0.00258  NA    
##  6      50        1  0.621    0.680    0.0590  3.73e-1 -0.0719   0.190    NA    
##  7     100        0  0.00175  0.0115   0.00974 2.99e-1 -0.00870  0.0282   NA    
##  8     100        1  0.586    0.556   -0.0302  6.40e-1 -0.158    0.0976   NA    
##  9     150        0  0.00403  0.0214   0.0174  2.48e-2  0.00221  0.0325   -0.183
## 10     150        1  0.579    0.525   -0.0538  4.09e-1 -0.183    0.0750   NA    
## 11     200        0  0.00558  0.00680  0.00122 8.48e-1 -0.0112   0.0137   NA    
## 12     200        1  0.574    0.570   -0.00377 9.55e-1 -0.136    0.128    NA    
## # … with abbreviated variable names ¹​in_effect, ²​mn_peak_cond1, ³​mn_peak_cond2
peak_df_1 |> 
  filter(in_effect == 1) |> 
  ggplot(aes(n_trial, group = n_trial))+
  geom_point(aes(y = mn_peak_cond1), colour = line_colours[1], shape = 15)+
  geom_point(aes(y = mn_peak_cond2), colour = line_colours[2], shape = 15)+
  geom_hline(aes(yintercept = 0), size = 1.05)+
  ggtitle("Mean Amplitude of ERP Peaks in two Conditions")+ 
  labs(subtitle = glue::glue("Conditions:, 
                            <b><span style='color:{line_colours[1]};'>Control Group </span></b>and 
                            <b><span style='color:{line_colours[2]};'>Alternative</span></b>"))+
  scale_y_continuous("Mean Amplitude (mV)")+
  scale_x_continuous("Sample Size")

  # geom_point(aes(y = mn_diff), colour = "grey", shape = 15)
  # geom_errorbar(aes(ymin = mn_peak_cond2 - ci_l, ymax = mn_diff + ci_u), colour = line_colours[2])

# ggsave(filename = here::here("images", "ch4_mean_peak_1.png"), width = 10, height = 8, dpi = 360)

peak_df_1 |> 
  filter(in_effect == 1) |> 
  ggplot(aes(n_trial, group = n_trial))+
  geom_point(aes(y = mn_diff), colour = line_colours[1], shape = 15)+
  geom_errorbar(aes(ymin = mn_diff - abs(ci_l), ymax = mn_diff + abs(ci_u)), colour = line_colours[1], width = 5)+
  # geom_point(aes(y = mn_peak_cond2), colour = line_colours[2], shape = 15)+
  geom_hline(aes(yintercept = 0), size = 1.05)+
  ggtitle("Mean Difference in Amplitude of ERP Peaks in Various Sample Sizes")+ 
  scale_y_continuous("Mean Amplitude (mV)")+
  scale_x_continuous("Sample Size")+
  labs(subtitle = "Error Bars Represent 95% Confidence Intervals of Mean Difference")

# ggsave(filename = here::here("images", "ch4_mean_peak_2.png"), width = 10, height = 8, dpi = 360)